function[u2]=propRS(u1,L,lambda,z)
% propagation - transfer function approach
% assumes same x and y side lengths and
% uniform sampling
% u1 - source plane field
% L- source and observation plane side length
% lambda - wavelength
% z- propagation distance
% u2 - observation plane field

[M,N]=size(u1);                         % get input field array size
dx=L/M;                                 % sample interval
k=2*pi/lambda;                          % wavenumber
x=-L/2:dx:L/2-dx;
y=x;
r=sqrt(z^2+x.^2+y.^2);

fx=-1/(2*dx):1/L:1/(2*dx)-1/L;          % freq coords
fy=fx;
[FX,FY]=meshgrid(fx,fx);

H=exp(1i*k*z*sqrt(1-(lambda.*FX).^2-(lambda.*FY).^2));   % transfer function
H=fftshift(H);                          % Shift transfer function
U1=fft2(fftshift(u1));                  % shift, fft src field

U2=U1.*H;                               % multiply
u2=ifftshift(ifft2(U2));                % inverse fft, center observation field
end